Lattice density functional theory at finite temperature with strongly 
density-dependent exchange-correlation potentials 
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The derivative discontinuity of the exchange-correlation (xc) energy at integer particle number is 
a property of the exact, unknown xc functional of density functional theory (DFT) which is absent 
in many popular local and semilocal approximations. In lattice DFT, approximations exist which 
exhibit a discontinuity in the xc potential at half filling. However, due to convergence problems of the 
Kohn-Sham (KS) self-consistency cycle, the use of these functionals is mostly restricted to situations 
where the local density is away from half filling. Here a numerical scheme for the self-consistent 
solution of the lattice KS Hamiltonian with a local xc potential with rapid (or quasi-discontinuous) 
density dependence is suggested. The problem is formulated in terms of finite-temperature DFT 
where the discontinuity in the xc potential emerges naturally in the limit of zero temperature. A 
simple parametrization is suggested for the xc potential of the uniform ID Hubbard model at finite 
temperature which is obtained from the solution of the thermodynamic Bethe ansatz. The feasibility 
of the numerical scheme is demonstrated by application to a model of fermionic atoms in a harmonic 
trap. The corresponding density profile exhibits a plateau of integer occupation at low temperatures 
which melts away for higher temperatures. 

PACS numbers: 71.15Mb,71.10Fd,71.10Pm,71.27+a 
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I. INTRODUCTION 

Originally, static (ground-state) density functional the- 
ory (DFT) has been formulated^ for many-electron sys- 
tems in the continuous space of three spatial dimensions 
with the electrons interacting via the Coulomb interac- 
tion. On the other hand, many phenomena in many- 
particle physics are studied in terms of model systems 
on discrete lattices, typically of tight-binding form. The 
one-dimensional Hubbard model or the Anderson impu- 
rity model are just the most prominent examples. Typi- 
cally, these models are studied with techniques different 
from DFT. However, DFT can be a useful tool for the in- 
vestigation of these models as well, especially when one 
wants to take into account the effects of non-uniform ex- 
ternal potentials^—. For example, cold atoms in optical 
lattices confined in an harmonic trap may very well be 
modelled by a lattice model with confining external po- 
tential where the particles interact through a Hubbard- 
type interaction!. 

The idea of formulating DFT for electrons on a lat- 
tice of sites has been pioneered by Schonhammer and 
coworkers^. As with usual DFT, the applicability and 
success of lattice DFT hinges on the availability of ap- 
proximations to the unknown exchange-correlation (xc) 
functional. Capelle and coworkers proposed a local func- 
tional for the xc energy per site based on the Bcthc- 
ansatz solution of the uniform ID Hubbard model at zero 
temperature^. Thus, for one-dimensional lattice mod- 
els, the ID Hubbard model takes the role of the uniform 
electron gas in continuum formulation of DFT as an ex- 
actly solvable model system which provides the essential 



input for the construction of the local approximation. 
In the same spirit, xc functionals have been suggested 
based on other Hubbard lattice models such as the two- 
dimensional hexagonal latticed or the simple cubic lat- 
tice in 3D±±. 

An interesting property of the Bethe-ansatz local den- 
sity approximation (BALDA) in ID (and also its coun- 
terparts for 2D or 3D lattice models) is its discontin- 
uous form of the xc potential at half filling or integer 
occupation^. Physically, this discontinuity is a direct 
consequence of the Mott-Hubbard gap of the Hubbard 
model while in the DFT context it is nothing but the 
well-known derivative discontinuity of the xc energy at 
integer particle number at zero temperature^. 

The BALDA has been successfully applied to spa- 
tially inhomogencous Hubbard supperlattices^, to cold 
fermionic atoms in a harmonic trap, both with repulsive^ 
and attractive^ electronic interaction and to the study of 
the static and dynamic linear density respons o 15 ' 16 . Ex- 
tensions of the BALDA have been suggested to systems 
in static magnetic fields^ and, in the adiabatic form, to 
the domain of time-dependent DFT— where it has been 
used to study the dynamics of finite Hubbard clusters. 
Recently the adiabatic BALDA has been applied to de- 
scribe the time evolution of trapped ID lattice fermions 
in the Mott insulator regime^. A modified version of the 
BALDA has been used in the study of time-dependent 
transport through an Anderson impuritj*2£ where the dis- 
continuity has been related to Coulomb blockade. 

From a physical point of view the discontinuity at inte- 
ger particle number is certainly a desirable property for 
an approximation to have, at least at zero temperature. 
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As we will argue below, the zero-temperature disconti- 
nuity may be viewed as the zero temperature limit of 
a continuous xc potential at finite temperature. From a 
practical point of view, if a local discontinuous (or rapidly 
varying) xc potential is used, one often faces convergence 
problems of the Kohn-Sham (KS) self-consistency cycled 
essentially whenever the local density is close to integer 
occupation. 

In the present work we propose a practical solution to 
this convergence problem by viewing it as an equivalent 
problem of finding the solution to a coupled set of non- 
linear equations. In Sec. [Til we start with a general dis- 
cussion of the KS self-consistency cycle and possible con- 
vergence problems when using KS potentials which vary 
rapidly for small variations in the density. The problem is 
illustrated explicitly on the simple, exactly soluble model 
system of a single interacting site in contact with a heat 
and particle bath. In Sec. IIII1 we then introduce the ID 
lattice models studied throughout this work and briefly 
summarize the idea of a local approximation for lattice 
models which has been discussed in the literature. We 
will work in the framework of finite-temperature DFT—, 
and Sec. [IV] is devoted to the construction of an approx- 
imate xc potential for this framework. We construct the 
xc potential of the uniform ID Hubbard model for fi- 
nite temperatures based on the thermodynamic Bethe 
ansatz. We provide a simple parametrization of this po- 
tential using insights gained from the simple single site 
model discussed earlier. In Sec. [V] we introduce our al- 
gorithm for the practical solution of the self-consistency 
problem which is based on a multi-dimensional bisection 
method. In Sec. IVII we show a numerical application 
of the method to the problem of interacting particles 
in a harmonic trap before we present our conclusions in 



Sec. lVIII In the Appendix we provide explicit expressions 
for the xc free energy per site for the simple parametriza- 
tion of the thermodynamic Bethe ansatz solutions to the 
uniform Hubbard model. 



II. KOHN-SHAM PROBLEM WITH RAPIDLY 
VARYING DENSITY FUNCTIONALS 

The implementation of DFT via the KS method gained 
enormous popularity because it reduces calculations of 
the density n(r) in a complicated strongly interacting 
system to computing n(r) for a reference system of non- 
interacting KS particles. The KS particles move in the 
presence of an effective potential v KS = v+vuxc [ n ]> where 
v is an external potential and i>hxc[?i] is the Hartree- 
cxchangc-corrclation (Hxc) potential which depends on 
the density and is adjusted sclf-consistently to reproduce 
the physical density distribution of the interacting sys- 
tem. The self-consistent nature of the KS problem makes 
it nonlinear and thus not absolutely trivial. In fact, the 
whole point of the present paper is to identify one of the 
potentially dangerous physical situations and to propose 
a recipe for its solution. 



A. KS self-consistency as a fixed point problem: 
the issue of convergence 

Assuming that the potential i>hxc[?i] as functional of 
the density is known, the general KS problem can be 
formulated as follows. We have to find a set of KS or- 
bitals <^ Q ) and KS energies e a by solving a one-particle 
stationary Schrodingcr equation 



(i + v + VHxc[n])v? (a) = e a tp 



(1) 



where t is the one-particle kinetic energy operator. As 
the operator in Eq. (fl]) depends on the density we need an 
additional "self-consistency equation" that relates the set 
of {^ a \ e a } to n(r). Obviously this equation is simply 
the standard definition of the density of noninteracting 
particles 



n(r) = 2£/(e Q ) 



p (Q) (r) 



(2) 



where the factor two comes from spin, /(w) = (1 + 
exp(/3(w — /i))) _1 is the Fermi distribution, f3 = 1/T is 
the inverse temperature, and \x is the chemical potential 
which is cither given externally or determined by fixing 
the total number of particles. 

Calculation of the density from Eqs. (HJ)-© is equiv- 
alent to finding a fixed point of a certain density func- 
tional. Indeed, the eigenvalue problem of Eq. (fT]) defines 
a map n h-> {tp^ a \e a } from the density to the set of KS 
eigenfunctions and eigenvalues, i. e., it determines the 
functionals and e a [n] . Inserting these functionals 

into Eq. ^ we obtain the following form 



2 £/(<>«[»]) 



G[r. 



(3) 



which is a typical fixed point problem for the functional 
G[n] on the right hand side. 

In practice, the self-consistent KS problem of Eqs. JT])- 
@, or equivalently the fixed point problem of Eq. (J3|), 
is commonly solved iteratively. In the simplest case one 
starts with some initial guess for the density and 
constructs a sequence of iterations as follows 



G[i 



,(o)l 



,(*) 



G[i 



,( fc -Dl 



(4) 



The limiting point of this sequence presumably gives a 
self-consistent solution of the KS equation 



n = lim n 

fc- 



(fe) 



(5) 



Unfortunately the assumed convergence cannot be guar- 
anteed in general, in spite of the fact that the original 
KS problem definitely has a unique solution. From the 
Banach fixed point theorem (the contraction mapping 
principle) we know that the sequence of Eq. Q does 
necessarily converge to a unique fixed point if the func- 
tional F[n] is contractive, i. e., if the following condition 
is satisfied 

||G[n] -G[n'] || <A||n-n'||, < A < 1, (6) 
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where || ... || means a properly chosen norm in the space 
of densities. Apparently this condition requires F[n] to 
be a sufficiently smooth functional of the density, which is 
not always the case. Moreover there are important phys- 
ical situations where the inequality of Eq. © is always 
violated. To understand this more clearly we estimate 
the left hand side of Eq. ([6]) for a small density variation 
n' = n + Sn with Sn <C n 



|G[n]-G[n']|| 



X— - — (n-n) (7) 
on 



where x is the density response function. Obviously 
the right hand side of Eq. ([7]) cannot be smaller than 
A 1 1 77, — n'\\ with < A < 1 if the Hxc potential is a 
rapidly varying functional of n, i. e., if Sv g^ c is large at 
least for some directions in density space. Physically this 
should always happen in systems composed of weakly 
coupled fragments if the number of particles in at least 
one of the fragments is close to an integer value. Then 
a density transfer to/from this fragment causes a strong 
variation of the potential. The origin of this behavior is in 
the famous discontinuity of the exact xc potential at inte- 
ger number of particles^. The most prominent examples 
of systems demonstrating such a behavior are molecules 
close to dissociation or strongly correlated solids near the 
Mott-Hubbard transition. In all those systems where the 
physics is governed by a nearly discontinuous xc poten- 
tial the standard iterative procedure of solving the KS 
equations will not converge. 

In the next subsection we explicitly illustrate the above 
general argument by considering a very simple model sys- 
tem - a single lattice site which can host at most two 
spin-1/2 fermions. The purpose for studying this model 
is twofold. Firstly, this is probably the only case where 
the exact xc potential can be found analytically for any 
temperature. The corresponding KS problem possesses 
an analytic solution and, because of its simple structure, 
clearly shows when and why the existing unique fixed 
point cannot be reached iteratively. Secondly, a single 
site DFT serves as a paradigmatic example for more gen- 
eral interacting lattice models. In fact, the analytic form 
of the single site Hxc potential will later be used to con- 
struct a simple parametrization for the Hxc potential of 
the uniform Hubbard model at finite temperatures. 



B. KS-DFT for a single site model 

Let us consider one single-orbital site in contact with 
a heat and particle bath at inverse temperature /? and 
chemical potential / i 22 i 23 . 

The Hamiltonian for this single site model (SSM) in 
the presence of an on-site interaction is given by 



Hssm = v n + U rio.tncj. 



(8) 



where vq is the on-site energy and U is the charging en- 
ergy, no ;(T and ho = J2a=t l ^o,<t are the operators for 



the on-site density with spin a and for the total density, 
respectively Similarly, for the non-interacting case the 
single-site Hamiltonian reads 



SSM — VsTlO 



(9) 



with on-site energy v s . The complete Fock space of both 
Hamiltonians is spanned by the states |0), | f), | \), and 
| fj,) with particle occupation of zero, one, and two. 
These states are both eigenstates of Hssm with eigen- 
values 0, vq, vq, and 2vq + U, as well as eigenstates of 
Hi SM with eigenvalues 0, v s , v s , and 2u s , respectively. 
For the single site model, the particle number operator is 
equal to the density operator, N = ho, and the density 
no = Tr {pho} for the interacting case then reads 



n 



2exp(-/3(t. - n)) + 2exp(-/?(2(t; - p) + U)) 

_£SSM 



(10) 



where 

z ssm = 1 + 2 exp (_ /3 ( Vo _ M )) + eX p(-/3(2(u - aO + CO) 

(11) 

is the grand canonical partition function. Eq. (|10[) only 
depends on the quantity vq = vq — [i and the function 
n o(^o) can be inverted explicitly leading to 



„ , s „ 1, Sn+ JSn 2 + e-P u {l ~ 8n 2 )\ 

Mno) = - U -j3 hl [ —n ) 

(12) 

where Sn = no — 1. 

Following the same lines, for the non-interacting case 
the density reads 

, _ 2cxp(-(3(v s - fi)) + 2exp(-/3(2(t; s - 

n - ^|SM V 1 "^ 

with the non-interacting partition function 

Z SSM = x + 2cx p(-/3(z; s - n)) + exp(-0(2(v, - /i))) • 

(14) 



Again, the density only depends on the quantity v s 
v s — /i and one can invert ng(u s ) to yield 



(15) 



with Sn s = rig — 1. 



The exact Hxc potential for the SSM can now easily 
be calculated by requiring that the interacting density 
equals the non- interacting one no = =: n and taking 
the difference of the two expressions (fTSji and (|T2"j) . i.e., 



4 s x f ( n> u, T) = v s (n) - v (n) = ^ + g(n - 1) (16) 



where 



9W = | + I to £±V^M (17) 



n 

FIG. 1. Hartree-exchange-correlation potential of the single- 
site model for different temperatures T = 1//3. Energies given 
in units of U. 

which is easily shown to be an odd function of its argu- 
ment, g{-x) = -g(x). 

In Fig. Q] we show the SSM Hxc potential v| s x ^(n) as 
function of the density for different temperatures. At 
low temperatures, Wy^(n) becomes an extremely rapidly 
varying function of n in the vicinity of n = 1 , approaching 
a step function with a step of height U at n = 1 in the 
limit of zero temperature. 

Now, having at hand the exact Hxc potential, we can 
study the KS problem. In the single site DFT, the general 
fixed point equation reduces to the following algebraic 
transcendental equation 

n = 2f(v + vl™(n)) = G(n). (18) 

This equation can be easily solved analytically. By con- 
struction, the solution to Eq. (|18[) simply returns the 
function n(T,n) defined by Eqs. {TDJ) and (TT]). In Fig.H 
we show the left (straight line) and right hand sides of 
Eq. (f!8|) for two different values of vq — \x- Obviously, for 
a given vq there is only one intersection between n and 
G(n), which means that the function G(n) always has 
only one fixed point. However, if the expected solution 
lies in the region of fast variation of i>hxc, be. n ~ 1, 
the fixed point cannot be reached by iterations. Inde- 
pendently of the choice of the initial guess, after a few 
iterations we enter a limiting cycle with the density end- 
lessly jumping between n w and n w 2. This behavior 
is generic for low enough temperatures, T <C U , and the 
chemical potential in the region vq < /i < vq + U, which 
are the conditions ensuring that vn xc (n) has a step-like 
form (see Fig. [TJ) and the physical on-site occupation is 
close to unity. Examples of the iterative cycle are indi- 
cated in Fig. [2] showing the convergence of the cycle of 
Eq. (j4|) for vo — fx = —1.25, and the lack of convergence 
for Vq — /U = —0.25. 

From the discussion in Sec. Ill Al it is clear that the 
same type of non-convergence of the KS iterative se- 
quence should occur in any system with a discontinu- 
ous/rapidly varying xc potential. In the rest of this paper 



QU 1 I 1 1 1 1 , 

0.5 1 1.5 2 

n 

FIG. 2. Left-hand side (dotted black line ) and r.h.s. of 
Eq. (|18|l for different values of vq—^l and temperature T = 0.1. 
(Energies given in units of U.) For given vq — /x, the self- 
consistent solution is given by the intersection of the corre- 
sponding G(n) with the straight line. For —U < vo — /i < 0, 
the iterative scheme n' fc+1 ' = G(n' fe '), indicated by the dash- 
dotted lines, does not converge. 



we study and solve this problem for lattice models where 
the discontinuity of vrxc reflects Mott-Hubbard correla- 
tions and can easily be captured at the level of a local 
density approximation. 



III. LATTICE DENSITY FUNCTIONAL 
THEORY 

A. Lattice DFT: Formalism and Model 

As a particular example of a lattice model, we consider 
one-dimensional, interacting many-electron systems on a 
tight-binding lattice described by the Hamiltonian 

N B -1 

H = -t j2 E (4<^+b- + H - c ) 

i=l <T=t4 

N s N„ 

i=l i—1 

where c\ a (cj )Cr ) are the fermionic creation (annihila- 
tion) operators for an electron with spin a at site i, 
hi >rT = cj a Ci. a and n, = $^o-=t j, fii,a are the operators 
for the density of electrons with spin a and for the total 
electron density at site i, respectively. The nearest neigh- 
bor hopping element is t, U is the Hubbard interaction. 
Vi is the external potential at site i and N s is the total 
number of sites. For simplicity, we consider systems in 
the absence of magnetic fields. For the grand-canonical 
ensemble, when the system is in contact with a heat bath 
at inverse temperature /? and a particle bath at chemical 
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potential [i, the statistical operator is 

p=^exp(-p(H-fiN)) (20) 

where N = 1S the operator for the total number 

of particles and the grand canonical partition function 

is Z = Tr |exp(— /3(H — fj,N))\ with the trace over all 
states of Fock space. In thermal equilibrium, an observ- 
able described by the operator A then takes the value 

A = Tr {M}- 

In the spirit of DFT at finite temperatures^, the 
Hamiltonian (|19[) is mapped onto the following Hamil- 
tonian of non-interacting electrons 

N s N„ 

i=l CT=f4 »=i 

where the effective single particle KS potential v^ s at 
site i is chosen such that the equilibrium density n, = 
Tr{p?ii} of the interacting Hamiltonian (fTT)| and the KS 
Hamiltonian (f2"Tj) are the same for all sites. The KS po- 
tential at site i then has the form 

v? s =v i+ v?™ (22) 

where m is the external potential at site i and, similarly, 
u^ xc is the Hxc potential at site i. In general, the Hxc 
potential at site i depends on the equilibrium density at 
all other sites, i.e., vf xc = vf xc ({rij}). Typically, how- 
ever, the exact form of the Hxc potential is unknown and 
one has to resort to approximations. Once an approxima- 
tion to v^ xc has been specified, the equilibrium density of 
the KS Hamiltonian H KS can be found by self-consistent 
solution of the KS equation on the lattice 

E(-^+«f%)^ a) =^ Q) (23) 
j'=i 

(with tij = t for j = i± 1 and ty = otherwise) together 
with 

n t = 2j2f(e a )\^ a) \ 2 ■ (24) 

a 

B. Local density approximations in the lattice DFT 

In the local approximation, the Hxc potential at site i 
only depends on the density at the same site i, v i xc ' oc = 
v Hxc( n i)- The functional dependence of Dg° c d (n) on the 
density is extracted from some interacting model system 
for which the exact solution can be constructed by ana- 
lytical and/or numerical techniques. Probably the most 
prominent example of such a functional for lattice-DFT is 
the local density approximation (LDA) based on Bcthe- 
ansatz solution of the uniform Hubbard model in ID 
(Bcthe-ansatz LDA, BALDA) at zero temperature^&2i. 



Strictly speaking, at zero temperature and exactly at 
half-filling (n = 1), the BALDA xc potential is not de- 
fined since the xc energy per particle is not diffcrcntiablc 
at this point. One pragmatic way around this mathe- 
matical problem is to smoothen the discontinuity in some 
ad-hoc manne r 19 i 20 . Alternatively one can construct xc 
functionals for finite temperature which approach a dis- 
continuous function in the zero temperature limit. We 
have already followed this route in Sec. lIIBl to formulate 
a single site DFT, and will pursue it further in Sec. IIVI 
for the ID Hubbard model. 

Although one can avoid the use of truly discontinuous 
xc potentials in this way, the resulting KS potentials will 
still be very rapidly varying functions of the density. This 
is exactly the property leading to a non-convergence of a 
simple iterative procedure, a fact which has been recog- 
nized in attempts to use the BALDA xc potential within 
the usual KS self-consistency cycled. 

IV. LOCAL APPROXIMATIONS AT FINITE 
TEMPERATURE 

In the present Section we propose several versions of a 
local functional at finite temperature for which the corre- 
sponding Hxc potentials exhibit rapid variations as func- 
tion of the density. This functional is based on the ther- 
modynamic Bethe ansatz (TBA) solution of the uniform 
Hubbard model in one dimension 2 ^ and is thus an exten- 
sion of the corresponding work at zero temperature^. 

We have numerically solved the coupled integral equa- 
tions of the TBA following Refs. I26K28I . For given inverse 
temperature /3, the density is calculated as a function of 
the chemical potential which can be inverted to give the 
chemical potential as function of the density. For the 
interacting and non-interacting cases these inverse func- 
tions are denoted as fJ-(n) and /J. s (n), respectively. From 
these two functions we obtain the density-dependent Hxc 
potential of the TBA as 

«™ A (n) = - /x s (n) (25) 

In Fig. [3] we show the density dependence of the TBA 
Hxc potentials for various values of the interaction U and 
various temperatures T = 1//3. At low temperatures and 
for sufficiently large values of U, the TBA Hxc potential 
V™c ( n ) exhibits rapid variations around half filling (n = 
1) as function of density. In the zero-temperature limit 
this feature reduces to a step whose height is given by 
the Mott-Hubbard gap. This gap can be expressed in 
terms of the parameters of the model as (from now on all 
energies are given in units of the hopping matrix element 
t unless otherwise noted) 

16 f°° \/x 2 - 1 

A o(U) = 77 / dx -^k 777T (26) 

y ' U Ji sinh(27ra/[f) v ; 

which is nothing but the derivative discontinuity of the 
uniform Hubbard model at half fillin g 12 i 24 . 
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n 



FIG. 3. Comparison of the fully numerical TBA Hxc po- 
tential with pOTBA parametrization of Eq. (|27[) for different 
interactions and two typical temperatures. 

Away from n = 1, even at low temperatures the 
dependence of u™ A (ri) on the density is rather slow 
and smooth. For high temperatures, the sharp feature 
around half filling is washed out. As a consequence of 
particle-hole symmetry, the TBA Hxc potential takes 
the value U/2 at n = 1 for all temperatures and ex- 
hibits a point symmetry around this point as function 
of density, i.e., v^(n) = U/2 + g TBA (n - 1) with 
<? TBA (— x) — —g TBA (x). Furthermore the values of the 
TBA Hxc potential at the endpoints of the density inter- 
val are Whxc A (°) = and v Hx C A ( 2 ) = u for au tempera- 
tures. 

We use these observations to design a hierarchy of an- 
alytic parametrizations of the fully numerical TBA Hxc 
potential which can easily be used in practical calcula- 
tions. In the construction of these parametrization we 
make use of the simple analytic form of the Hxc poten- 
tial of the single site model discussed in Sec. IIIBI 



A. Lowest level single-site-motivated 
parametrization of the numerical TBA: pOTBA 

Our simplest functional is aimed at reproducing the 
main qualitative features of the full numerical TBA based 
Hxc potential. These are the point symmetry of the 
function w™c A ( n ) (reflecting the electron- hole symme- 
try), and the step structure at n = 1, which gradually 
washes out at higher temperatures. 

Precisely this pattern is also observed in the Hxc po- 
tential of the single-site model: in the zero-temperature 
limit, the Hxc potential u™(n) has a step of height 



U. This almost discontinuous feature at low tempera- 
tures crosses over to a smooth one at high temperatures. 
Hence we will use the analytic form of to mimic the 
step feature in our "lowest level" parametrization of the 
TBA Hxc potential. We adopt the simplest possible way 
to reproduce the correct low temperature amplitude of 
the step. Namely, in the function v^^(n,U,T), defined 
in Eqs. (|16|) - (fTT)) . the parameter U will be replaced by the 
zero-temperature Mott- Hubbard gap Aq(U) of Eq. (|26p . 
The reduction of the gap automatically reduces the value 
of the potential at n — 2 from the exact value of U down 
to Ao(U). This unwanted behavior is corrected by adding 
a proper linear function that also ensures the right point 
symmetry of the Hxc potential. Putting all these argu- 
ments together we propose the following simple zero-level 
parametrization for the TBA Hxc potential (pOTBA) 

<™ A W = U -f {U) n + v s ™(n,A (U),T). (27) 

The Hxc potential defined by this equation is shown in 
Fig. [3] together with the full numerical «™ c A (n). Wc 
clearly sec that for all T and U the parametrization pro- 
posed in Eq. (|27|) overall agrees reasonably well with 
u™ A (ri). The maximal deviations never exceed unity 
(i.e., the value of the hopping integral t). Obviously the 
simple linear form of the first term of Eq. (f2"Tf is not flex- 
ible enough to reproduce the detailed features of the full 
TBA Hxc potential away from the step (see Fig. [3J . How- 
ever, as long as the difference between the parametriza- 
tion and the full TBA Hxc potential are small compared 
to t, these inaccuracies are, in most practical cases, of 
little consequence for the solution of the KS equations, 
as will be confirmed in Sec. I VII It is also worth not- 
ing that a reasonably accurate practical approximation 
of Eq. (f2"T|) does not actually require the solution of the 
TBA equation. The only input we used was a zero- 
temperature Mott-Hubbard gap and general symmetry 
arguments. This observation can be useful to construct 
local approximations for more complicated, e. g. multi- 
dimensional, lattice models for which no exact solutions 
are available. 

Apparently there are cases when the fine structure of 
the density distribution cannot fully be captured within 
our simplest zero-level parametrization pOTBA defined 
by Eq. If!??)) . Therefore it is desirable to design a re- 
fined parametrization which further reduces the devia- 
tion from the numerical TBA potential. Two successive 
refinements of the "first-level" (plTBA), and of the "sec- 
ond level" (p2TBA) are described in the next two sub- 
sections. 



B. First-level refined parametrization correcting 
the temperature dependence: plTBA 

Fig. [3] clearly shows that the temperature dependence 
of our simple pOTBA potential Eq. ((27)) is not per- 
fect. At higher temperatures the step in the function 
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w Hxc BA ( n ) washes out too fast as compared to the nu- 
merical w™ A (n). There is an obvious physical reason 
for this deficiency. When the temperature T increases 
and becomes larger than unity (in units of the hopping 
integral t), the kinetic energy contribution to the parti- 
tion function becomes less and less important. Therefore 
at T > 1, and independently of U, the system should 
behave more or less like a collection of independent sites 
with the Hxc potential given by the pure SSM expression 

of Eqs. (HSD-dIZJ>- 

In our first-level refinement (plTBA) we take into ac- 
count this physics by replacing Aq(U) in Eq. (j2"T)l with a 
"temperature-dependent gap" Ai({7, T) 



v 



PlTBA^ _ 
Hx 



(n) 



U-Ax(U,T) 



+ v*? K f(n,A 1 (U,T),T). 



(28) 

The function &i(U,T) reduces to A (U) at T < 1 and 
approaches U in the opposite limit of T > 1. The two 
limits are connected by a smooth function which is de- 
termined by comparison with the numerical TBA data. 
We have found that the following Padc-likc form does the 
required job 



a m rp\ A Q (U)+a^(T)T + UT 2 
Ai(t7,JJ= x^f 2 



(29) 



with 



(T) 



'T + a 



(i) 



T 2 + 



,(1) 



Here = 0.95, 



-0.08, 4 X) = 0.13. 



From Fig. @] we see that the first-level paramctrization 
plTBA, Eqs. P8" ]) -(|29" ]) . produces an Hxc potential which 
is practically indistinguishable from the full numerical 
f™ c A (n), provided that T is not too small. However, 
there are still some deviations in the low-temperature 
regime. This point is addressed at the last step in our 
three-level hierarchy of parametrizations. 



C. Second-level refinement — the best analytic fit 
to numerical TBA: p2TBA 

The only feature missing in the first-level plTBA 
paramctrization is a low-temperature nonlinearity of 
u™ A ()j) away from the half filling (see upper panel of 
Fig.@|. Physically the nonlinearity should be attributed 
to a nontrivial density of states in the Hubbard bands. 
According to our experience the remaining inaccuracy in 
the Hxc potential has practically no effect on the den- 
sity distribution. On the other hand, we cannot exclude 
that in some situations (steep external potentials and / or 
small number of particles) the low-temperature inaccu- 
racy of ^Hxc BA ( n ) w ^ produce visible (thought definitely 
not large) errors in the density. To avoid such problems 
we go to the last step in our hierarchy and introduce a 
nonlinear correction term. A very satisfactory fit to the 
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FIG. 4. Comparison of the fully numerical TBA Hxc potential 
with the first-level refined plTBA parametrization of Eq. 
for two typical temperatures. 



numerical v]^^(n) can be achieved with the following 
(p2TBA) form 



p2TBA i \ _ 



7 Hxc W - 2 " ^ UHxc 

-A(U,T)sm[2ir(n- 1)] - B(U,T) sin[7r(n- 1)].(30) 

We note that the analytic form of the correction term 
in Eq. (|3Tj)) automatically preserves the point symmetry 
of the potential and the exact values at the end points, 
n = and n = 2. Note also that we use a new function 
A 2 (Z7,T) defined by 

A (U)+aW(T)T + UT 2 



U-A 2 {U,T) 



+ <f(n,A 2 (C/,T),T) 



A 2 (U,T) 



1 + T 2 



(31) 



with 



,(2) 



,( 2 



(T) 



T 



(2) 



T 



(2) 



,(2) 



-0.28, a 2 2) 



2.2, and 



and the coefficients arc a\ 

(2) 

a 3 = 0.50. As before, in the zero-temperature limit 
A2 (t/,T) reduces to the correct zero- temperature gap 
Aq(JJ) while in the high-tempcrature limit it becomes 
U. The interaction and temperature dependent coeffi- 
cients A(U,T) and B(U,T) in Eq. ([30) are parametrized 
as follows 



A(U,T) 



AijTp 2 
U 2 + A 2 (T) 



with 



Ai(T) 



^11 
T' 2 + A 12 



B(U,T) 



A 2 (T) 



BjjTp 2 
U 2 + B 2 (T) 



A 2 



T 2 + A 22 



,(32) 



(33) 
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n 



FIG. 5. Comparison of the fully numerical TBA Hxc poten- 
tial with the second-level refined p2TBA parametrization of 
Eq. (|30[) for two typical temperatures. 

BUT) = — , B 2 (T) = — ^ — . (34) 

1 ' T 2 + B 12 ' 21 ' T 2 + B 22 V ' 

Here A n = 0.09, A 12 = 0.25, A 21 = 705.5, A 22 = 24.95 
and B n = 0.05, B 12 = 0.13, B 21 = 0.65, B 22 = 0.01. 
The accuracy of this parametrization can be appreciated 
in Fig. [5J at the scale of the plot, the parametrization 
p2TBA and the full numerical TBA Hxc potentials are 
essentially indistinguishable. 

Our parametrization of the finite-temperature TBA 
results generalizes earlier parametrization a 4,30 valid for 
zero temperature. Close comparison of the Hxc poten- 
tial of our p2TBA parametrization of Eq. (|3T)|) in the 
zero-temperature limit with the FVC parametrization of 
Ref. [10 and exact zero-temperature results reveals that 
the p2TBA parametrization in some density ranges can 
be marginally less accurate than FVC. As pointed out 
before, however, in contrast to FVC our parametrization 
by construction incorporates the exact zero-temperature 
gap. Since the most prominent feature of the Hxc poten- 
tial as function of density is precisely the discontinuity 
at half- filling, i.e., the zero-temperature gap, in some sit- 
uations the incorrect gap of FVC might lead to spurious 
features as will be shown below. 



V. SELF-CONSISTENCY WITH RAPIDLY 
VARYING FUNCTIONALS USING BISECTION 

In the previous Section we presented a hierarchy of ex- 
plicit local approximations for the Hxc potential of lat- 
tice DFT, which at low temperatures are rapidly varying 
functions of the density close to half- filling. In Sec. iHlwc 



have pointed out the difficulties in converging the usual 
self-consistency cycle for solving the KS equation ([1) with 
such functionals. In the present Section we show how to 
avoid the convergence problem and present a numerically 
feasible algorithm to obtain the self-consistent solution 
based on bisection techniques. 

We begin by writing again the self-consistency equa- 
tion for the density at site i, Eq. (|24|) . in a form of a 
fixed point problem, i. c., making explicit its dependence 
on the densities at all other sites: 

n i =J2f(^(n))\ V ( r ) (n)\ 2 = 0^] (35) 

a 

where 

n = (ni,n 2 ,...,njvj (36) 

and the orbitals are calculated from the KS equation ((23]) 
using the KS potential 

v? s (n i )=v i +vg£(n i ). (37) 

The set of equations ([33)) for i 6 1, . . . , N s constitute a 
coupled set of N s nonlinear equations for the N s ground- 
state densities m, i = 1, ...,N S . As we have argued 
in Sec. |TTJ the numerical solution of these equations by 
plain iterations produces a non-converging sequence. Ob- 
viously more elaborate iterative schemes, such as the 
Newton-Rhapson method, which incorporates informa- 
tion on derivatives of the equations with respect to 
the unknown variables are also not appropriate because 
those derivatives may become very large (in the low- 
temperature regime). This again leads to convergence 
problems in the iterative solution of the coupled nonlin- 
ear equations. 

Here we propose a solution scheme based on bisection. 
The main idea of our algorithm is inspired by the exactly 
solvable single-site KS problem described in Sec. IIIBI 
Therefore we first explain it for this simple, but very 
illuminating case. Instead of iteratively searching for the 
fixed point of function G(n) in Eq. (|18p we rewrite this 
equation as 

n - G(n) = 0, (38) 

and search for zeros of the left hand side. From the 
uniqueness of the solution we know that there is only 
one zero, and, by construction, the l.h.s. of Eq. ((38]) 
has different signs at the end points, n = and n . = 2, 
of the density interval. Therefore the standard bisection 
method^ is applicable. Starting from the end points and 
using bisections we can bracket the solution to any de- 
sired accuracy. 

For the general lattice KS problem we need to solve a 
system of Eqs. ([3"5|) . In this case the following straight- 
forward multidimensional generalization of the standard 
bisection method can be used. We start with an initial 
guess n- ^ for the densities at sites i 6 {2,..., N s }. Wc 
then define the density vector 

nW = (n 1 ,n 2 0) ,...,n<$l) , (39) 
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FIG. 6. Solution lines of rii — Gi(ni,ri2) = 0, i = 1,2, for 
the symmetric three-site problem. The intersection of these 
two lines is the solution to the system of Eqs. (|42[) - (|43l) . The 
parameters are U = 8, fj. = 0.5, Vi = 2, vi = 0, and T = 0.5. 



insert this density vector in the rhs of Eq. (|55|) for i = 1 
and solve the resulting nonlinear equation for the density 



ni 



with the usual ID bisection method. Then we 



choose the density vector 



n 



(2) = 



,712,71-3 



(0) 



(0) 



insert it into Eq. (|3"5j) for i = 2 and solve for 712 
again by bisection. In the next step we take 



(3) - L« „w 



n 2 ,n 3 ,n 4 



(0) 



(0) 



(40) 



(41) 



and solve Eq. (|35[) for i = 3 for 713 = We continue 
the procedure until we have exhausted the N s equations 
(f3"5"1) . Then we start the cycle all over but now with the 
initial guess for the density n^ Ns+1 ^ as Eq. (|39[) but with 

the nf 3 ^ replaced by nj . The whole process is continued 
until convergence is achieved. 

In Fig. [6] we illustrate the iterative procedure for a 
symmetric three-site problem, i.e., for V\ = W3 and by 
symmetry also v^ s = v^ s . We then need to solve two 
coupled nonlinear equations of the form 



n-i - Gi(m, n 2 ) = 0, 
n 2 - G 2 (ni, n 2 ) = 0. 



(42) 
(43) 



In Fig. [S] we show the {n\, n 2 ) plane and zero lines for the 
l.h.s. of Eqs. g2|) and g3|). The overall solution of the 
problem is given by the intersection of these two lines. 
Furthermore, the dashed line in Fig. [5] shows how our 
iterative scheme converges to this solution. 

From this example we can also understand that in some 
cases the proposed scheme, as is, may not converge. If the 
zero lines close to the intersection become just straight 



lines then, depending on the slopes of these lines, the it- 
erative scheme may follow a rectangular path encircling 
the solution point but never reaching it. However, in 
this case one additional Newton-Rhapson step (using in- 
formation on the derivatives) will directly lead to the 
solution point. 

As is common in DFT, we have used the densities as 
fundamental variables in Eq. (|3"5|) . It is also possible to 
implement the bisection scheme in terms of the KS po- 
tentials. To this end we write the vector of Hxc potentials 
at sites i as 



VHx 



XC, 1 J * * * 3 ^Hxc,iV s ) 



(44) 



and consider both the KS orbitals and eigenvalues as 
functions of this vector, i.e., = v?^( v Hxc) and 

e (a) — e^^vhxc)- Therefore, also the density at site 
i can be considered a function of vjjxc through 



ni(vHxc) = X^( £Q ( VHxc ))l'^ a) ( VHxc )l 2 



(45) 



For our local approximations to the Hxc potentials the 
set of nonlinear equations to be solved by bisection then 
becomes 



^x°c d K(v H xc)) 



(46) 



with the local density ti^vhxc) given by Eq. (|45[) . In or- 
der to solve this set of equations we proceed in an anal- 
ogous way to the one described above for the density- 
based procedure but now treating the local Hxc poten- 
tials VHxc.i at site i as the unknowns to be determined. 



VI. NUMERICAL APPLICATIONS 

In this Section we present some numerical exam- 
ples demonstrating successful applications of our self- 
consistency algorithm as well as the accuracy of the local 
approximation using our different parametrizations. 

We illustrate our theoretical developments by calcu- 
lating the density distribution of particles confined by an 
external harmonic potential of the form 



Vi = V cxt (i-ioY 



(47) 



where V cx t is the strength of the trapping potential and 
io = (N s + l)/2 (we take V ext in units of the hopping 
parameter t). The Hubbard chain with a superimposed 
harmonic potential is commonly used to model the be- 
havior of cold fcrmionic gases in ID optical lattices— 
Therefore our results below have a clear relevance for the 
physics of cold trapped atoms. However, for our present 
illustrative purposes, the choice of this particular system 
is related to one of its specific features, namely the pos- 
sible coexistence of the Mott insulator phase around the 
center of the trap and the metallic phase at the trap's 
edges^ i 31 ' 32 . In the Mott phase the density is pinned at 
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FIG. 7. Density distribution of N — 70 particles in the 
harmonic potential at T — (upper panel: U = 2, lower 
panel: U = 8). The strength of the external potential is 
Kxt = 2.5 x l(T :i . We compare the exact DMRG results 
against DFT results based on different choices for uhxc- They 
are the parametrization of Ref. [33 (labeled as FVC), our 
lowest-level parametrization pOTBA of Eq. (|27|l (labeled as 
pOTBA), and the full numerical TBA potential (labeled as 
TBA), respectively. The results obtained with the p2TBA 
potential of Eq. p0[) are indistinguishable from the full TBA 
results. 



n = 1 which shows up as an extended plateau in the den- 
sity distribution. Therefore in trapped systems the ap- 
pearance of the Mott insulator phase becomes detectable 
within the "density-only" DFT concept. On the other 
hand, at the level of DFT functionals the Mott physics 
is solely related to the discontinuity of the xc potential. 
Thus the Hubbard model with a harmonic confinement 
is perfectly suited for demonstrating the working power 
of our algorithm as well as the performance of different 
paramctrizations for the Hxc potentials. 

As a first example we study a system with N = 70 
particles on N s = 100 sites in the presence of the poten- 
tial given by Eq. g7]) with V ext = 2.5 x 10~ 3 . In Fig. 
we show self-consistent densities for two different values 
of the Hubbard interaction, U — 2 and U — 8, evalu- 
ated at zero temperature. Except for the DMRG results 
which denotes numerically exact reference results from 
density matrix renormalization group calculations^—, 
all other calculations result from self-consistent DFT cal- 
culations with local approximations to the Hxc poten- 
tial. FVC denotes the zero-temperature BALDA using 
the parametrization of Ref. [3(1 pOTBA denotes results 
obtained with our low-level parametrization of Eq. (|27p , 
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FIG. 8. Melting of the Mott phase at finite temperature. The 
calculations are based on the p2TBA parametrization which 
produces the results identical to those of the full numerical 
TBA. The parameters are (7 = 8, Kxt = 2.5 x 10" 3 , N = 70. 



the results obtained from our second-level parametriza- 
tion of Eq. ((30)) arc indistinguishable from those using 
the exact numerical TBA as input in the local approxi- 
mation. 

We see that the density profiles are quite similar 
in the different approaches. For U = 2 the FVC 
parametrization exhibits two spurious density plateaus 
around i ~ 30 and i k 70 which are due to the fact 
that this parametrization does not incorporate the zero- 
temperature derivative discontinuity exactly but only ap- 
proximately. We also see some small differences between 
the pOTBA and the TBA results in the flanks of the den- 
sity profile. For stronger interaction, U = 8, the density 
exhibits an extended plateau of value unity over a wide 
range of sites in the center of the well which physically 
corresponds to the local, incompressible Mott phase. On 
the other hand, in the density functional picture this 
plateau is a direct consequence of the extremely rapid 
variation of the Hxc potentials as function of the density. 
The small difference in the pOTBA and TBA densities 
is related to the nonlinearity of Hxc potential away from 
the half filling (see Sec. lIVC) . This deficiency is corrected 
in our second level parametrization p2TB A of Eq. (|30[) . 
As a result the density calculated with the p2TBA Hxc 
potential is completely indistinguishable from that ob- 
tained using the full numerical TBA potential. This also 
holds true for all interactions, temperatures and trapping 
potentials we have tried. Hence in practice in all figures 
the results denoted as TBA have been actually produced 
using the p2TBA potential defined after Eq. ([501) • At 
strictly zero temperature, the KS potential has a real 
discontinuity at integer filling n = 1 and, therefore, it is 
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FIG. 9. Comparison of the density distributions at T = 
obtained using Hartree-Fock approximation (HFA), Thomas- 
Fermi approximation (TFA) against DFT results based on 
the TBA potential (TBA). The parameters are the same as 
in Fig. [7J 




FIG. 10. Density distribution of TV = 10 particles in a har- 
monic trap at temperature T — in Hartree-Fock (HFA) and 
Thomas- Fermi approximations (TFA), compared to DFT re- 
sults based on TBA and DMRG. Upper panel: U — 2, lower 
panel: (7 = 8. The external potential is V cx t = 2.5 x 10 -2 . 



undefined in this point. For practically solving Eqs. (|35[) - 
([57]) in this case we replace the discontinuity by a linear 
function when n G [1 — An, 1 + An] with An = 10 -3 . 
Further decreasing An will not hamper the convergency 
and the final results. In general, the densities from the 
DFT calculations are remarkably close to the numerically 
exact quantum Monte-car lo/DMRG result&^i^. 

Unlike earlier work^ 3 ^, our parametrization is valid 
for arbitrary temperature and thus allows to study fi- 
nite temperature effects. As a first application of this 
feature, in Fig. [5] we show how the density plateau due 
to the local Mott phase "melts away" when increasing 
the temperature. 

We have also calculated (Fig. ^ the density pro- 
files at zero temperature in the Hartree-Fock (HFA) 
and Thomas- Fermi approximation (TFA). Here, by TFA 
we mean that the non-interacting kinetic energy is not 
treated exactly as in the Kohn-Sham scheme but at the 
level of a local approximation. It is important to note, 
however, that exchange-correlation effects are also in- 
cluded at the level of the local density approximation (in 
Ref. i this approximation has been denoted as "total- 
energy LDA (TLDA)"). The top panel of Fig. M shows 
that for small values of U both the HFA and the TFA 
give a reasonably accurate density profile when compared 
to DMRG. For larger values of U (U = 8, lower panel of 
Fig. however, the situation is different: HFA com- 
pletely misses the development of the density plateau 
while TFA (due to the local Hxc potential) does ex- 



hibit this plateau. Moreover, in HFA the density is more 
spread out as compared to the exact result. 

The differences between TFA and both the DMRG 
as well as the full KS results with the TBA functional 
are more pronounced for smaller number of particles. In 
Fig. [TO] we show the density profile for N = 10 electrons 
on N s = 40 sites in the harmonic external potential of 
Eq. (gTJ) with V cxt = 2.5 x 10~ 2 . Not surprisingly, the 
TFA approximation completely misses the quantum os- 
cillations in the density profile both for small and large 
values of U. In contrast, HFA captures these oscillations 
well for small U since the kinetic energy is treated at a 
quantum level. In contrast, for U = 8 the correlation ef- 
fects are too strong to be captured by HFA which, again, 
shows a density distribution which is too spread out. On 
the other hand, the KS calculation using the TBA again 
is in extremely good agreement with the DMRG results. 
The lower panel of Fig. [TU] also shows that for this case 
our simplest pOTBA parametrization of the TBA results 
can sometimes lead to inaccuracies. Finally, in Fig. 111! 
solving the KS equations with the TBA functional, we 
show the effect of temperature on the density profile: one 
can see that the quantum oscillations due to the relatively 
small number of particles seen for zero temperature are 
quickly suppressed when increasing the temperature. 



12 




FIG. 11. Effects of finite temperatures on the density profile 
of N = 10 particles in a harmonic trap. The calculation is 
based on the full numerical TBA. The parameters are U = 8, 
Kxt = 2.5 x 1(T 2 . 



VII. CONCLUSIONS 



In this work we have suggested a way to deal with 
convergence problems in self-consistent KS calculations 
when dealing with (local) approximations for the xc 
potential which exhibit rapid variations as function of 
the density. Working in the framework of lattice DFT, 
we have formulated the KS self-consistency cycle as a 
fixed point problem and shown that for rapidly vary- 
ing functionals the fixed point in the usual prcedure is 
not attractive. Instead, we rephrased the search for 
the self-consistent KS potential in terms of finding the 
roots of a set of coupled nonlinear equations. To find 
these solutions, we then suggested an iterative algorithm 
based on successive application of the well-known bi- 
section method for finding roots of nonlinear equations 
in one dimension. The scheme has been successfully 
tested for model systems of electrons in a harmonic trap 
interacting via a Hubbard interaction. We have used 
a newly designed local approximation for the xc func- 
tional based on the thermodynamic Bcthc ansatz solu- 
tion of the uniform Hubbard model. Based on these 
results we constructed simple, yet accurate parametriza- 
tions for arbitrary temperatures, thus generalizing earlier 
parametrizations suggested for zero temperature. This 
paves the way for further investigations on the perfor- 
mance of finite-temperature DFT for one-dimensional 
lattice models. 
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Appendix: Hartree-exchange-correlation free energy 
per site for the TBA parametrizations 



In Sec. |IV]we have suggested a local density approxi- 
mation for one-dimensional lattice systems based on the 
numerical solution of the TBA for the uniform Hubbard 
model. The parametrizations we proposed were con- 
structed from insights gained on the numerical results 
for the Hxc potential. However, typically the construc- 
tion of local DFT approximations starts by modelling the 
xc energy per site and the corresponding xc potential is 
then obtained by differentiation. In this Appendix we 
derive the expressions for the Hxc free energies per site 
for the different parametrizations suggested in Sec. IIVI 

We start with the derivation of the exact LDA Hxc 
free energy per site at finite temperature, expressed in 
terms of general thermodynamic quantities. As usual, a 
crucial ingredient of the general construction of LDA is 
a reference system of interacting particles with uniform 
density for which the grand canonical potential per site, 
is written as function of the chemical potential fi. 
The partition function for the reference system is 

Z(p) = exp(-0fiO*)) (A.l) 

from which we can derive the density as function of //, 

an(fi) idiaz{n) 



dfj, 



(A.2) 



This function can be inverted to give the chemical po- 
tential as function of density, /j, = n(n). By Legcndrc 
transformation we can then obtain the free energy per 
site of the reference system as function of density as 



F(n) = fi(/i(n)) + n(n)n . 



(A.3) 



We can repeat the same steps for the corresponding non- 
interacting reference with the same uniform density n 
with grand canonical potential per site Q s ((j, s ) and the 
corresponding expression for the free energy per site 

F s (n) = Q s O s (n)) + n a (n)n . ( A.4) 

The Hxc free energy per site is then simply given by 

F Hxc {n)=F{n)~F s (n) 

1 f Z(fi s (n) +v Hxc (n))\ 
= { ZMn)) ) +nmM (A ' 5) 
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where we have defined the Hxc potential as 

UHxc(n) = fj,(n) - n a (n) (A. 6) 

(see also Eq. ([25])). 

Applying Eq. (|A.5[) to the single-site model of Section 
III BL we obtain for the Hxc free energy per site of that 
model 



F^(n,U,T) = nv^(n,U,T) 

1 ( Z^(^{n,T) + v^(n,U,T)y 

{ Z SSM (/1 SSM (n|T) ) 



where the Hxc potential of the SSM model, 
Vx^f(n,U,T)), is given by Eq. and we have 

made explicit the dependence on the temperature T 
and the on-site interaction U. It remains to find the 
dependence of Z SSM and Zf SM on the density. This 
can be done by performing the program described 
above with the interacting and non-interacting partition 
functions of the SSM model given in Eqs. (Til"]) and (fT4"| . 
Actually, we have already calculated the dependence of 
the chemical potential on the density (see Eq. ([T5])) 



fif M (n,T) = -v.(n) = ^ In 



/3 V 2 - n 



(A. 



which, when inserted back into Eq. (TH| . yields the non- 
interacting partition function of the SSM model in terms 
of the density 



2 - n 



(A.9) 



Finally, the interacting partition function of the SSM 
model (see Eq. (JT])) may be written in terms of the den- 
sity as 



Z SSM (n, U, T) = 1 + JlL exp (K^K ^ T)) 
■exp(/3(2z;f x f(n,[/,T)-[/)) (A.10) 



(2-n) 2 

which leads to the final result for the Hxc free energy per 



site of the SSM model 



SSM , 



Hxc (n,U,T) = nv^(n,U,T) 



In 



1 - 



-^)exp(/34 s x f(n,t/,T)) 



-^exp(/3(24 s x f(n,C/,T)-[/)) 



(A.11) 



By construction, the derivative of this expression with 
respect to the density yields the SSM Hxc potential of 
Eq. dUl). 

Using this result it is now easy to express the Hxc 
free energies per site for our different parametrizations 
of the TBA results. For the lowest-level parametrization 
(pOTBA) of the TBA the resulting expression reads 



F, 



pOTBA 



(n, U, T) = U A / U) n* + F^(n, A (U), T) 

(A.12) 

where Ao(u) is the exact zero-temperature gap of the 
uniform Hubbard model given by Eq. (|2l)]) . For the first- 
level refinement (plTBA), it has the same form except 
that A (U) is replaced by Ai({/,T) of Eq. (gU), i.e., 



plTBA 



Hxc 



(n,U,T) = 



U-Ax(U,T) 



+F^ c M (n,A 1 ([/,T),T) .(A.13) 



Finally, for the second refined parametrization (p2TBA) 
we have 



F 



^ BA (n,U,T) = 



+F^(n,A 2 (U,T),T) + 
B(U,T) 



U-A 2 (U,T) 



A(U,T) 
2tt 



(cos(?r(n - 1)) + 1) 



(cos(2tt(?i- 1)) - 1) 
(A.14) 



with A2(U,T) given by Eq. (|3"Tj) and the functions 
A(U,T) and B(U,T) given by Eq. (32). 

It is worth mentioning that for zero temperature, since 
A 2 (U,0) = Ai(U,0) = A (U), the contributions to the 
Hxc free energy per site coming from the single-site model 
in all our parametrizations (|A.12[) - (|A.14[) reduce to the 
same limit 



? SSM 
Hxc 



(n, A (U), 0) = A (U)(n - l)6(n - 1) (A.15) 



where 0(x) is the Heaviside step function. 
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